Blantyre Air Quality Modelling

Load required libraries

library(tidyverse)

── Attaching core tidyverse packages ───────────────────────────────────────────────────────────────────────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.2     ✔ tibble    3.2.1
✔ lubridate 1.9.4     ✔ tidyr     1.3.1
✔ purrr     1.0.4     
── Conflicts ─────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors

library(here)

G3;here() starts at /Users/peter.macpherson/Library/Mobile Documents/com~apple~CloudDocs/Documents/Documents - Peter’s Mac mini/Projects/AQ-MTB_2025-05-14/Work/Data/aq_mtb
g

library(brms)

G3;Loading required package: Rcpp
gG3;Loading 'brms' package (version 2.22.0). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').
gG3;
Attaching package: ‘brms’

gG3;The following object is masked from ‘package:stats’:

    ar

g

library(tidybayes)

G3;
Attaching package: ‘tidybayes’

gG3;The following objects are masked from ‘package:brms’:

    dstudent_t, pstudent_t, qstudent_t, rstudent_t

g

library(sf)

G3;Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE
g

library(mgcv)

G3;Loading required package: nlme
gG3;
Attaching package: ‘nlme’

gG3;The following object is masked from ‘package:dplyr’:

    collapse

gG3;This is mgcv 1.9-3. For overview type 'help("mgcv-package")'.
gG3;
Attaching package: ‘mgcv’

gG3;The following objects are masked from ‘package:brms’:

    s, t2

g

library(priorsense)
library(osmdata)

G3;Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright
g

library(nngeo) 

G3;Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
g

Import the modelling dataset - cleaned and prepped in the 2025-05-12_aq_clean.Rmd script

aq_model_data <- read_rds("input_data/aq_model_data.rds")

#also load the scale clusters
load("input_data/scale_72_clusters.rda") #SCALE clusters

#and the grid data
read_rds("input_data/mw_100m_cropped.rds")

$mwi_ppp_2020.tif
               [,1]         [,2]         [,3]         [,4]         [,5]         [,6]         [,7]         [,8]         [,9]        [,10]
  [1,]  0.045165006  0.046947237  0.055178735  0.047102101  0.043815471  0.041929662  0.039621390   0.03764743 3.824615e-02 3.529998e-02
  [2,]  0.054517046  0.047781128  0.052100033  0.046856463  0.042057369  0.028544975  0.027005022   0.02849834 2.889967e-02 2.749119e-02
  [3,]  0.045086842  0.047777772  0.041951492  0.034255959  0.029440563  0.024227845  0.016665278   0.02430716 7.012391e-03 1.374129e-02
  [4,]  0.035714723  0.039766397  0.035008151  0.032815386  0.030017985  0.028784296  0.029581040   0.02444118 1.111266e-02 4.598498e-03
  [5,]  0.015322668  0.018739052  0.018032806  0.018089430  0.017559310  0.018649841  0.017364351   0.01674390 1.384361e-02 5.014876e-03
              [,11]        [,12]        [,13]        [,14]        [,15]        [,16]        [,17]        [,18]        [,19]        [,20]
  [1,] 3.581402e-02  0.033903446  0.028697260  0.025856392 9.365267e-03 4.969447e-03 4.880766e-03 6.018387e-03 7.693967e-03 1.687009e-02
  [2,] 3.080107e-02  0.031454172  0.036637034  0.028348375 1.773694e-02 5.785842e-03 6.604420e-03 6.460932e-03 8.471724e-03 1.088048e-02
  [3,] 2.850383e-02  0.033486295  0.035977926  0.021789778 7.277420e-03 5.582869e-03 6.420635e-03 6.296951e-03 8.607163e-03 7.565869e-03
  [4,] 4.251544e-03  0.014496404  0.024309691  0.025809046 5.967221e-03 5.781442e-03 5.793448e-03 8.620093e-03 8.066839e-03 7.239441e-03
  [5,] 3.449674e-03  0.004229320  0.027064180  0.021735243 2.350955e-02 6.558061e-03 5.777395e-03 6.241026e-03 6.318555e-03 7.906212e-03
              [,21]        [,22]        [,23]        [,24]        [,25]        [,26]        [,27]        [,28]        [,29]        [,30]
  [1,] 7.483703e-03 1.134669e-02 2.477654e-02 1.747209e-02   0.02301022   0.03323989   0.03443765   0.04373282   0.03239537   0.03268054
  [2,] 6.492021e-03 1.798505e-02 3.820698e-02 3.622469e-02   0.04873810   0.06264146   0.04502079   0.04902079   0.03930074   0.04743655
  [3,] 1.413537e-02 2.868031e-02 4.856283e-02 7.635828e-02   0.10526711   0.10480442   0.09942765   0.08063150   0.07026784   0.07916132
  [4,] 3.765117e-02 4.892825e-02 8.083706e-02 6.265849e-02   0.06352694   0.10991651   0.11036292   0.10957718   0.09102185   0.09063222
  [5,] 4.181869e-02 6.396030e-02 6.208089e-02 5.934627e-02   0.06683331   0.11021938   0.11196631   0.10069560   0.08459925   0.09193873
              [,31]        [,32]        [,33]        [,34]        [,35]        [,36]        [,37]        [,38]        [,39]        [,40]
  [1,]   0.01216105 6.328823e-03   0.01512775   0.02006114 3.101150e-02   0.05473484   0.06160709   0.07647948   0.07574075   0.07371373
  [2,]   0.04709647 1.577023e-02   0.01547693   0.01512168 9.752575e-03   0.03868737   0.05658274   0.06831920   0.07236640   0.06479138
  [3,]   0.08027354 4.902577e-02   0.04347426   0.04361450 2.506110e-02   0.02168380   0.03363772   0.07003722   0.10318663   0.07466165
  [4,]   0.10364193 1.078961e-01   0.06888389   0.06897457 7.109120e-02   0.06687172   0.04877206   0.06761945   0.06558794   0.05342549
  [5,]   0.10197861 1.167887e-01   0.11340023   0.07957450 8.267088e-02   0.07813363   0.05316366   0.03048283   0.04607106   0.05050959
             [,41]       [,42]       [,43]        [,44]        [,45]        [,46]        [,47]        [,48]        [,49]        [,50]
  [1,]  0.09937067  0.09119225  0.08600094   0.09103850   0.09306126   0.17105776   0.16657399   0.19017223   0.18780756   0.25865054
  [2,]  0.08613592  0.09276834  0.08768278   0.09463739   0.09931904   0.20142318   0.34862801   0.24239205   0.20016080   0.14733949
  [3,]  0.08796223  0.08529528  0.05861057   0.08760586   0.17597923   0.18829565   0.30547777   0.18160681   0.10488693   0.05029128
  [4,]  0.07356653  0.04595139  0.06095412   0.07659898   0.20957065   0.22002560   0.26929030   0.13957511   0.13037296   0.12920222
  [5,]  0.05924074  0.12837498  0.10924210   0.04965207   0.03692738   0.13063972   0.13985389   0.22550066   0.22139461   0.23821150
              [,51]        [,52]        [,53]        [,54]        [,55]        [,56]        [,57]        [,58]        [,59]        [,60]
  [1,]   0.18521543 1.718710e-01   0.17261122 1.592420e-01   0.13821790 2.070713e-01 2.180091e-01 2.510735e-01 3.347855e-01   0.42194489
  [2,]   0.25576282 1.704141e-01   0.16087675 6.019334e-02   0.05911760 2.086819e-01 2.241761e-01 2.393655e-01 2.442953e-01   0.28544879
  [3,]   0.08434317 1.395781e-01   0.03054703 2.303165e-02   0.02988590 1.481759e-01 1.958029e-01 2.119675e-01 2.289245e-01   0.22243160
  [4,]   0.12479196 1.226818e-01   0.12564214 1.318173e-01   0.13020913 1.327822e-01 2.432093e-01 2.560078e-01 2.778168e-01   0.30719724
  [5,]   0.22799887 1.374678e-01   0.14205216 1.374896e-01   0.14270458 2.208541e-01 2.599400e-01 2.764551e-01 4.585261e-01   0.61238015
              [,61]       [,62]        [,63]       [,64]        [,65]       [,66]       [,67]       [,68]       [,69]       [,70]       [,71]
  [1,]   0.20092252   0.2094545   0.23813984   0.2600536   0.28496569   0.3515317   0.3495567   3.0187855   3.4292254   4.3702903  10.2430563
  [2,]   0.37174022   0.2000832   0.20835581   0.2420900   0.25386712   0.3162111   0.3505953   2.8476219   2.6364815   3.5290518   4.9171147
  [3,]   0.31947485   0.3825941   0.37987089   0.2292890   0.24251479   0.2557049   2.4554725   2.3886611   1.5429624   2.8309567   3.7996442
  [4,]   0.34511879   0.5266556   0.32725015   0.2879633   0.32016167   0.8491008   2.2822711   3.1561983   4.8729448   4.3276281   5.2610259
  [5,]   1.26680827   1.3391285   1.04976392   0.3362455   0.29695830   0.2514963   0.7431583   6.0843654   6.2483072   5.5046158   5.5970216
             [,72]       [,73]       [,74]       [,75]       [,76]       [,77]       [,78]       [,79]       [,80]       [,81]       [,82]
  [1,]  12.3174086  20.4170475   5.9770074   5.8101201   6.1680431   7.0288363   2.9920485   2.7822750   3.2919111   2.7028601   2.7344890
  [2,]  11.4855623   3.7609022   3.5631394   2.8713186   2.5883286   6.6772809   6.6999898   6.8617449   3.1268249   2.7026587   2.1398206
  [3,]   4.0249171   2.3710566   2.3969679   1.5280255   1.3292254   1.3552140   6.3890018   6.8920774   5.9485264   6.3546324   2.1996837
  [4,]   1.3237036   1.1526742   2.2051387   1.6958532   5.7821097   1.2228683   2.1854522   6.9267850   5.9227381   5.7966475   2.1904876
  [5,]   1.9726193   1.2746783   1.6196227   1.4951222   5.7002840   1.2784475   5.3949475   2.2073100   2.0879290   2.0115442   1.1903656
             [,83]       [,84]       [,85]       [,86]       [,87]       [,88]       [,89]       [,90]       [,91]       [,92]       [,93]
  [1,]   2.8262646   3.1302779   3.6975005   2.1694636   2.2199211   2.3952386   2.0133886   1.0360730   6.0406532   6.6095033   6.3419876
  [2,]   1.7706543   1.8591148   2.6041703   1.9700989   2.0043545   2.1224627  12.0202541  10.6813793   5.7113476   5.8507247   5.5529652
  [3,]   0.7009542   0.5428279   1.0550138   0.9149451   9.6139355  10.1034241   9.7381792   7.6451182   5.3572125   5.1382504   5.1880946
  [4,]   0.6239597   0.5024059   0.9726501   0.9038598   9.6681395   7.2640862   6.9079909   5.7228518   6.3007116   6.8418317   7.1582618
  [5,]   0.4369978   0.4899155   0.8078130   0.7381263   8.1139145   6.2780352   5.4850063   5.4929547   6.6182308   8.1480541   7.9830341
            [,94]      [,95]      [,96]      [,97]      [,98]      [,99]     [,100]     [,101]     [,102]     [,103]     [,104]      [,105]
  [1,]  6.0880580  4.9141769  5.6921430  5.6045480  5.6583200  5.7986665  4.8305206  4.5144963  4.3507361  4.8232894  4.1951013   4.6705604
  [2,]  5.3891821  5.2441683  4.3157148  4.7750669  4.9417787  5.7972479  5.7227149  6.0235286  5.1872697  5.3412123  5.1507735   4.7766061
  [3,]  5.2370167  5.1144838  4.4092488  4.2801752  4.8011594  5.8726711  5.7670879  6.0404143  6.7926741  6.6716657  5.2236753   5.0849514
  [4,]  7.2552085  6.9078765  5.8223286  5.2468567  5.3122015  6.6020727  6.6301870  6.3754029  7.2572618  7.1930847  7.2747545   7.0348992
  [5,]  7.6431541  7.6273088  6.5719180  5.9137011  5.9537210  6.5491629  6.6482563  6.1040545  6.8009992  7.0888548  7.2325339   6.9409246
            [,106]      [,107]      [,108]      [,109]      [,110]      [,111]      [,112]      [,113]      [,114]      [,115]      [,116]
  [1,]   4.0513868   4.5824718  31.0151443  35.9615288  31.9965153  33.5947914  34.8950462  33.6029587  32.1246338  31.0912933  13.1747513
  [2,]   4.3982000   4.5165482  30.4003181  31.6037254  30.7493134  31.7781143  35.9739647  37.8679657  37.2548790  15.6511335  14.6784821
  [3,]   5.1197200   4.4855256  33.6331902  39.4704590  25.6988716  33.0550232  13.5869970  13.1186905  13.0235996  13.0852079  16.8719463
  [4,]   7.5529518   6.8021359  38.6025314  38.8714218  16.4733143  13.1496153  11.0804739   9.1725912  11.2854414  11.4492359  14.8278837
  [5,]   7.0277257   7.7998891   6.9491796  14.6560345  17.5421658  16.1886940  11.8365650  11.4074469  13.3405104  13.5668907  14.1344957
            [,117]      [,118]      [,119]      [,120]      [,121]      [,122]      [,123]      [,124]      [,125]      [,126]      [,127]
  [1,]  17.4172173  14.5559740  10.0551395   0.9203696   0.4967222   0.4107275   0.3978702   0.4675624   0.7120231   0.5725706   0.5796408
  [2,]  16.7342072  11.8646288  11.6098852   0.8615871   0.6930147   0.6287227   0.4980069   0.7469794   0.5880303   0.4945861   0.4778275
  [3,]  17.6651497  16.5033340  14.2913799   1.1762685   1.0414619   0.9750723   0.9179240   0.6648465   0.6006990   0.4671351   0.4315493
  [4,]  15.0737171  14.5026913   1.4341162   1.7666528   1.2385327   1.2163831   1.1903278   1.0431690   0.6815187   0.6674115   0.6116236
  [5,]  14.0739212   0.9036093   1.0228646   1.7256023   1.1906258   1.1074116   1.0897969   0.8469090   0.8302435   0.7722517   0.7019941
            [,128]      [,129]      [,130]     [,131]     [,132]     [,133]     [,134]     [,135]     [,136]     [,137]     [,138]     [,139]
  [1,]   0.5407016   5.6824121   5.5070438   5.393878   5.314744   5.378588   5.516953   5.472218   6.871101   7.120472   2.490248   1.913114
  [2,]   0.4782741   0.5772134   5.7760925   5.619320   5.583245   5.598059   5.671175   6.422924   7.114538   6.954295   2.615454   2.365523
  [3,]   0.4778309   6.0533495   5.8488727   5.827873   5.702469   5.692585   2.055359   2.058331   2.605898   2.534510   2.526598   1.925066
  [4,]   6.0364213   6.8703146   6.9667211   6.728249   2.507176   2.410894   2.355129   2.363863   2.335734   2.941173   2.280494   2.165629
  [5,]   7.4033751   6.2904220   6.0554729   5.708776   2.480082   2.428432   2.440901   2.502323   2.452162   2.456221   2.371370   2.334642
           [,140]     [,141]     [,142]     [,143]     [,144]     [,145]     [,146]     [,147]     [,148]     [,149]     [,150]      [,151]
  [1,]   1.890159   1.864571   1.824639   1.813932   1.855493   1.862693   1.844257   1.884519   1.856309   1.911597   1.900523   1.8835605
  [2,]   1.930770   1.902335   1.879635   1.835053   1.837441   1.848145   1.740682   1.844983   1.900461   1.933071   1.847667   1.8871386
  [3,]   1.854827   1.917260   1.922128   1.888102   1.913544   1.864133   1.771954   1.868454   1.890933   1.831353   1.933270   1.9211165
  [4,]   2.187655   2.135860   2.126067   2.083079   2.073362   2.087825   2.053761   2.042283   2.012455   2.036653   2.052104   1.9967043
  [5,]   2.811734   2.288969   2.239703   2.176805   2.207034   2.119645   2.070684   2.082571   1.989093   2.088601   2.078266   2.0321856
            [,152]      [,153]      [,154]      [,155]     [,156]     [,157]     [,158]     [,159]     [,160]      [,161]      [,162]     [,163]
  [1,]   1.8760191   4.2506528   1.8467858   1.9328438   4.556878   4.581502   5.842960   5.753704   4.339655   5.3159437   3.6693192  3.5465922
  [2,]   1.8791692   4.2897072   4.4196887   4.4844222   4.750791   5.543705   5.940609   5.749850   4.410793   5.4709144   5.1160870  3.6569257
  [3,]   1.8962737   4.3825073   4.1637969   4.6726365   4.731542   4.696297   5.871382   5.844340   4.788056   5.8306942   3.9437084  3.8018162
  [4,]   2.6364172   2.5760887   4.5244684   4.6207194   5.156001   5.343033   5.413106   5.551797   5.356392   3.1116621   3.1301842  3.1519930
  [5,]   2.6615610   2.5835884   2.7367139   2.7678025   5.207396   5.004921   4.730000   4.490225   6.635791   3.2056251   3.5913785  3.6644924
           [,164]     [,165]     [,166]     [,167]     [,168]     [,169]    [,170]    [,171]    [,172]    [,173]    [,174]    [,175]     [,176]
  [1,]  3.4214435  3.2192101  3.0704648  3.2135775  4.2359557  4.4121046  4.179620  5.511585  7.292460  4.900838  4.766211  4.537093  3.3284004
  [2,]  3.6246648  3.5184062  3.2872734  3.1027408  4.3657980  5.0197420  6.747860  8.048831  7.825340  7.773330  4.906329  4.761338  3.9066331
  [3,]  3.5798521  3.6366291  3.9564095  5.8185658  6.5702181  6.7610564  6.643492  7.909804  7.872230  8.015193  5.043985  5.217367  6.0537314
  [4,]  2.3774395  3.6038818  3.4744658  4.0604777  4.2596602  4.2969036  4.415190  5.865735  7.650286  7.599905  7.653657  5.146161  5.9680948
  [5,]  3.3508301  3.3644121  3.6195219  4.1825180  4.3101077  4.4454837  5.979053  9.620084 15.794232 10.375804  9.111597  5.311432  5.9213576
           [,177]
  [1,]  3.8432529
  [2,]  9.6162081
  [3,] 11.1199265
  [4,] 11.5795832
  [5,] 11.8466797
 [ reached 'max' / getOption("max.print") -- omitted 172 rows ]

attr(,"dimensions")
$x
$from
[1] 2755

$to
[1] 2931

$offset
[1] 32.67125

$delta
[1] 0.0008333333

$refsys
Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        MEMBER["World Geodetic System 1984 (G2296)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]

$point
[1] FALSE

$values
NULL

attr(,"class")
[1] "dimension"

$y
$from
[1] 7627

$to
[1] 7803

$offset
[1] -9.363749

$delta
[1] -0.0008333333

$refsys
Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        MEMBER["World Geodetic System 1984 (G2296)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]

$point
[1] FALSE

$values
NULL

attr(,"class")
[1] "dimension"

attr(,"raster")
$affine
[1] 0 0

$dimensions
[1] "x" "y"

$curvilinear
[1] FALSE

$blocksizes
        x y
[1,] 3892 1

attr(,"class")
[1] "stars_raster"
attr(,"class")
[1] "dimensions"
attr(,"class")
[1] "stars"

Model 1:

“We fitted a spatially smooth regression model for log-transformed PM2.5 concentrations using a Gaussian process (GP) smooth over spatial coordinates (x, y) to capture flexible spatial trends. Seasonal variation was modelled using a Fourier series expansion of day-of-year up to the 3rd harmonic (i.e., sine and cosine terms for annual, semi-annual, and tri-annual cycles) to account for complex seasonal patterns. The model was fitted in a Bayesian framework using brms with a Gaussian likelihood and the cmdstanr backend.”

Note, priors are still causing issues for me, so we will comment out now, and use the defauly stan weakly informative priors

Priors

priors <- c(
  prior(normal(3.4, 1), class = "Intercept"),
  prior(normal(0,1), class = "b"),
  prior(exponential(1), class = "sigma"),
  prior(exponential(1), class = "sdgp"),
  prior(normal(0, 1), class = "lscale", lb = 0)
)

Fit prior only model

Now model with data.

pp_check(m1)

G3;Using 10 posterior draws for ppc type 'dens_overlay' by default.
g

Compare prior-only model to model with data using priorsense package:

Predictions by space and time. Here, although we only have data within clusters, we will predict for cells outside of clusters based on covariates. We will also predict for the month for whcih we do not have data.

#set_up prediction date frame
nd_m1 <- st_as_sf(mw_100m_cropped)  %>%
  mutate(centroid = st_centroid(geometry)) %>%
  mutate(x = st_coordinates(centroid)[, 1],
         y = st_coordinates(centroid)[, 2]) %>%
  select(-centroid)  %>%
  st_drop_geometry() %>%
  rename(pop_density = 1) %>%
  mutate(pop_density_km2 = pop_density / 0.01)

G1;H1;Errorh: object 'mw_100m_cropped' not found
Error during wrapup: not that many frames on the stack
Error: no more error handlers available (recursive errors?); invoking 'abort' restart
g

Now just draw 50 random coordinates, and predict by time to see whether we captured the trend.

#sample coordinate points from the household dataset
#we will fix the small number sampled later!
set.seed(123) 
sampled_points_m1 <- aq_model_data %>%
  sample_n(50) %>%
  select(x, y)

#generate DOY 0–365 and calculate Fourier terms
doy_grid <- tibble(doy = 0:365) %>%
  mutate(
    sin_doy1 = sin(2 * pi * doy / 365),
    cos_doy1 = cos(2 * pi * doy / 365),
    sin_doy2 = sin(4 * pi * doy / 365),
    cos_doy2 = cos(4 * pi * doy / 365),
    sin_doy3 = sin(6 * pi * doy / 365),
    cos_doy3 = cos(6 * pi * doy / 365)
  )

#expand grid across sampled locations
prediction_df_m1 <- sampled_points_m1 %>%
  crossing(doy_grid)

#add predictions
preds_m1 <- add_epred_draws(object = m1, newdata = prediction_df_m1)

#summarise
preds_m1_sum <- preds_m1 %>%
  ungroup() %>%
  group_by(doy) %>%
  summarise(.epred_mean = mean(.epred),
            .lower = quantile(.epred, probs=0.025),
            .upper = quantile(.epred, probs = 0.975))

#GEt the observed data
observed_data <- aq_model_data %>%
  select(mean_doy, log_pm2_5) 

#plot
preds_m1_sum %>%
  ggplot() +
  geom_ribbon(aes(x=doy, ymin=.lower, ymax = .upper), fill="steelblue", alpha=0.3) +
  geom_line(aes(x=doy, y=.epred_mean)) +
  geom_jitter(data = observed_data, aes(x = mean_doy, y = log_pm2_5),
              color = "darkred", alpha = 0.6, width = 0.5, height = 0.0, size = 1.2) +
  labs(title = "Model-estimated log(PM2.5) with empirical measurements",
       subtitle = "Model predictions with 95% CrI and observed data points",
       x = "Day of year",
       y = "log(PM2.5)",
       caption = "Modelled estimates restricted to within clusters") +
  theme_ggdist() +
  theme(panel.background = element_rect(colour = "grey78"),
        legend.position = "none")

Model 2: With covariates

Now we include grid level covariates of distance to the road, population density, and building footprint percent, along witht the mean temp and huidity on the day of measurement (from the purple air monitors)

Again, priors to be fixed later

Predictions by space and time

Sample coordinates, and plot over time

#sample coordinate points
#will sort out the small number later
set.seed(123)
sampled_points_m2 <- aq_model_data %>%
  sample_n(50) %>%
  select(x, y, mean_temp_c, mean_current_humidity, pop_density_km2, building_coverage_pct, dist_to_road_m, grid_id, log_pm2_5)

#generate DOY 0–365 and calculate Fourier terms
doy_grid <- tibble(doy = seq(0,365, by=1)) %>%
  mutate(
    sin_doy1 = sin(2 * pi * doy / 365),
    cos_doy1 = cos(2 * pi * doy / 365),
    sin_doy2 = sin(4 * pi * doy / 365),
    cos_doy2 = cos(4 * pi * doy / 365),
    sin_doy3 = sin(6 * pi * doy / 365),
    cos_doy3 = cos(6 * pi * doy / 365)
  )

#expand grid across sampled locations
prediction_df_m2 <- sampled_points_m2 %>%
  crossing(doy_grid)

#add predictions
preds_m2 <- add_epred_draws(object = m2, newdata = prediction_df_m2, ndraws = 100)

#observed data for plotting
observed_data <- aq_model_data %>%
  select(mean_doy, log_pm2_5) %>%
  mutate(mean_doy = round(mean_doy))

#plot
preds_m2 %>%
  ungroup() %>%
  ggplot(aes(x=doy)) +
  stat_lineribbon(aes(y=.epred), .width=0.95) + 
  geom_jitter(data = observed_data, aes(x = mean_doy, y = log_pm2_5),
              color = "darkred", alpha = 0.6, width = 0.5, height = 0.0, size = 1.2) +
  scale_fill_brewer() +
  labs(title = "Model-estimated log(PM2.5) with empirical measurements",
       subtitle = "Model predictions with 95% CrI and observed data points",
       x = "Day of year",
       y = "log(PM2.5)",
       caption = "Modelled estimates restricted to within clusters") +
  theme_ggdist() +
  theme(panel.background = element_rect(colour = "grey78"),
        legend.position = "none")

Compare models, using LOO CV

library(loo)

loo_m1 <- loo(m1)
loo_m2 <- loo(m2)

loo_compare(loo_m1, loo_m2)

TODO

• OTHER OUTCOME MEASURES (PM1, PM10, NO)
• CAN WE LINK TO THE STATIC AQ MONITORS (OUTDOORS)
• PRIORS
• OTHER COVARIATES
• COMPARE MODELS WITH DIFFERENT BASIS FUNCTIONS FOR GP AND SPLINES
• EXCEEDANCES
• LINK TO MTB INFECTION PREVALENCE DATA, AND MODEL